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Abstract 



The Immersed Boundary method has evolved into one of the most useful computa- 
I't'I ■ tional methods in studying fluid structure interaction. On the other hand, the Immersed 

' Boundary method is also known to suffer from a severe timestep stability restriction when 

• . using an explicit time discretization. In this paper, we propose several efficient semi- 

^ ' implicit schemes to remove this stiffness from the Immersed Boundary method for the 

two-dimensional Stokes flow. First, we obtain a novel unconditionally stable semi-implicit 
discretization for the immersed boundary problem. Using this unconditionally stable dis- 
^ . cretization as a building block, we derive several efhcient semi-implicit schemes for the 

00 ' immersed boundary problem by applying the Small Scale Decomposition to this uncondi- 

0^ ■ tionally stable discretization. Our stability analysis and extensive numerical experiments 

' show that our semi-implicit schemes offer much better stability property than the explicit 

scheme. Unlike other implicit or semi-implicit schemes proposed in the literature, our semi- 
implicit schemes can be solved explicitly in the spectral space. Thus the computational 
. cost of our semi-implicit schemes is comparable to that of an explicit scheme, but with a 

' much better stability property. 

> '■ 

: 1 Introduction 

, The Immersed Boundary method was originally introduced by Peskin in the 1970's to model 

the flow around heart valves. Now it has evolved into a general useful method in studying the 
motion of one or more massless, elastic surface immersed in an incompressible, viscous fluid, 
particularly in biofluid dynamics problems where complex geometries and immersed elastic 
membranes are present. The method has been successfully applied to a variety of problems 
including blood flow in the heart [251 HI IIZl IISl ESI IHl [20] , vibrations of the cochlear basilar 
membrane [21 [8] , platelet aggregation during clotting [71 |3l] , aquatic locomotion [5l |6l jTH |35l [3] , 
flow with suspended particles [5l[3T], and inset flight [2H I22j. We refer to [27j for an extensive 
list of applications. 

The Immersed Boundary method employs a uniform Eulerian grid over the entire domain 
to describe the velocity field of the fluid and a Lagrangian description for the immersed elastic 
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structure. The force generated by the elastic structure drives the fluid and the fluid moves 
the elastic structure. This interaction is expressed in terms of the spreading and interpolation 
operations by use of smoothed Delta functions. 

One of the main difficulties that the Immersed Boundary method encounters is that it 
suffers from a severe timestep restriction in order to keep the stability |27t [5^ [50] . This has 
been the major limitation of the Immersed Boundary method. This restriction is typically 
much more severe than the one that would be imposed from using an explicit discretization for 
the convection term in the Navier-Stokes equation. The instability is known to arise from large 
boundary force and small viscosity [32]. Much effort has been made to remove this restriction. 
Some implicit and semi-implicit methods have been proposed in the literature [33l |23l [T5] . 
Despite of these efforts, the timestep restriction has not been resolved satisfactorily. The 
computational cost of using an implicit or semi-implicit scheme is still too high to be effective 
in a practical computation. To date, almost all practical computations using the immersed 
boundary method have been performed using an explicit discretization. 

In this paper, we develop several efficient semi-implicit schemes to compute the motion 
of an elastic interface immersed in a two-dimensional, incompressible Stokes flow. There are 
several important ingredients in deriving our semi-implicit schemes. The first one is to use the 
arclength and tangent angle formulation to describe the dynamics of the immersed interface 
[9]. We remark that Ceniceros and Roma have also used the arclength and tangent angle 
formulation to alleviate the stiffness of the viscous vortex sheet with surface tension in Q . The 
second one is to obtain an unconditionally stable semi-implicit discretization of the immersed 
boundary problem. Throughout this paper, we use the term "stability" to mean that the energy 
norm of the solution can be bounded in terms of the energy norm of the initial data, which is 
a weaker result than proving that the difference between two solutions in the energy norm can 
be bounded in terms of the energy norm of their difference at time zero. The third ingredient 
is to perform Small Scale Decomposition to the unconditionally stable discretization to obtain 
our efficient semi-implicit schemes. An important feature of our small scale decomposition 
is that the leading order term, which is to be discretized implicitly, can be expressed as a 
convolution operator. This property enables us to solve for the implicit solution explicitly 
using the Fourier transformation. Thus, the computational cost of our semi-implicit schemes 
is comparable to that of an explicit method. This offers a significant computational saving in 
using the Immersed Boundary method. 

The Small Scale Decomposition was first developed by Hou, Lowengrub and Shelley [9l[T0]. 
They applied this method to remove the stiffness from interfacial fiow with surface tension, 
which has proved to be very successful. Due to the coupling between the elastic boundary 
with the fluid, it is more difficult to remove the stiffness induced by the elastic force in the 
Immersed Boundary method. To remove the stiffness in the Immersed Boundary method, 
we need to decouple the stiffness induced by the elastic force from the fluid flow in such a 
way that the resulting semi-implicit discretization is still unconditionally stable. This is ac- 
complished by using a semi-implicit discretization which preserves certain important solution 
structures which exist at the continuous level. Without obtaining this unconditionally stable 
semi-implicit discretization, a straightforward application of the Small Scale Decomposition 
to the Immersed Boundary method would not provide an efficient semi-implicit scheme with 
the desirable stability property. Very recently, Newren et al. have obtained an uncondition- 
ally stable discretization for linear force in [23]. However, they did not perform Small Scale 
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Decomposition to their unconditionally stable discretization. As we will demonstrate in this 
paper, the unconditionally stable semi- implicit discretization without using the Small Scale 
Decomposition is still very expensive and the gain over the explicit discretization is quite 
limited. 

We develop several efficient semi-implicit schemes for both the steady Stokes flow and the 
unsteady Stokes flow respectively. In both cases, our semi-implicit schemes work very well. 
In the steady Stokes flow, we also develop a fourth order semi-implicit scheme by using the 
integral factor method. For the unsteady Stokes flow, we develop a second order semi-implicit 
method by combining our Small Scale Decomposition with a well known second order temporal 
discretization [131 127]- To illustrate the stability properties of our semi-implicit schemes, we 
apply our methods to several prototype problems and test our schemes for a range of elastic 
coefficients and viscosity coefficients. Our numerical results confirm that the semi-implicit 
schemes remove the high order stability constraint induced by the elastic force. In the case 
of unsteady Stokes equation, we also confirm the second order accuracy of our semi-implicit 
scheme. 

This paper is organized as follows. First, we review the classical formulation of the Im- 
mersed Boundary method in Section 2. Then, we introduce the arclength and tangent angle 
formulation in Section 3. In Section 4, we describe the spatial discretization of the Immersed 
Boundary method. In Section 5-6, we develop the numerical schemes for steady Stokes flow 
and unsteady Stokes flow respectively. The numerical results are presented in Section 7. Our 
numerical studies will focus on the stability restriction and computational cost of our methods. 
Some concluding remarks are given in Section 8. 



2 Review of the Immersed Boundary method 

For simplicity, we just consider a viscous incompressible fluid in a two dimensional domain 
0, containing an immersed massless elastic boundary in the form of a closed simple curve 
r. The configuration of the boundary is given in a parametric form: X(a, t),0 < a < L;,, 
X(0, = X(Lb,t), a tracks a material point of the boundary. We consider only the Stokes 
equations in this paper and would neglect the convection term. Then the governing equations 
are given as follows: 

P-Ql = -Vp + /xAu + f(x,t) , (1) 
V-u = 0, (2) 

(9X 

— (a,t) = u(X(a,t),t), (3) 

where u is the fluid velocity, p is the pressure, p and fj, are constant fluid density and viscosity, 
f(x, t) is the force density, which is not zero only on the boundary and which is inflnite there. 
The force density can be expressed as below 

f (x, t) = f^' F(q, t)6{x - X(a, t))da, (4) 
Jo 

5 denotes the two-dimensional Dirac delta function and 

F(a,t) = ^(rr), (5) 



3 



T 



ax 



da 



The choice of function T in this paper is computed by Hook's law 

9X 



Sb 



da 
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(6) 



(7) 



where Si, is the elastic coefficient of the boundary, and r is the unit tangent vector along the 
boundary, which is defined as 



ax 

ds , 



ax 



ds 



(8) 



This choice of force density has been used widely in the literature in both computational and 
theoretical studies [12] . [29] . [33] . 

We can rewrite ([3]) in the following way: 



ax, X 



/ u(x, t)5(x — X(a, t))(ix. 



(9) 



Next, we introduce the spreading and interpolation operations. The spreading and interpola- 
tion operators are defined as follows: 

L(X)(5(a))(x) = y ff(a)5(x-X(a,f))da, (10) 
L*(X)(u(x))(a) = / u(x)5(x-X(a,t))(ix . (11) 

It is easy to show that L and L* are adjoint operators: 

<u(x),L(X)(5(a)) >f, 
= y u(x) (7(a)(5(x — X(a, t))(ia^ dx 

= / / u(x)(5r(a)(^(x — X(q!, t))(iadx 

= u(-x.)g[a)6{:x. — X.{a,t))d:x.da 

Jr Jn 

= J g{a) (^J u{x)5{x — 'K{a,t))dx^ da 

= <L*(X)(n(x)),5(a) >r, 
where the inner product are defined as follows: 

<u,v>fi = / u(x)7;(x)dx, 
< f,g >r = / f{a)g{a)da. 



(12) 

(13) 
(14) 

Equations ([I]),© are the familiar Stokes equations of viscous incompressible fluid. Equa- 
tions ([3]),(I3]) represent the interaction of the fluid and the elastic boundary. The elastic bound- 
ary applies the force to the fluid, the fluid carries the immersed boundary, and the force density 
is determined by the configuration of the boundary. 
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3 The arclength-tangent angle formulation 



In studying the evolution of a curve, it is useful to represent the curve by its tangent angle 6 
and local arclength derivative s^. Previously, Hou, Lowengrub and Shelley [S] exploited this 
formulation and combined it with a so-called "Small Scale Decomposition" reformulation to 
remove the stiffness induced by surface tension. 

Consider the evolution of a simply closed curve T with known normal and tangent velocity 
fields, U,V. Assume the curve is represented by X(a,t),a G [OjL^]. We define the arclength 
derivative, Sa, and the tangent vector, 9, as follows 

Sa{a,t) = \Xa{a,t)\, (15) 
{xaio:,t),yaia,t)) = Sa{a,t){cos 9{a,t),sm9{a,t)). (16) 

The closed curve T evolves according to 

— = n{X,t) = Un + VT, (17) 

where r and n are the unit tangent and normal vectors of the curve respectively. According 
to the Frenet formula, we have ^ = kn, ^ = —kr, here s is the arclength variable. It is easy 
to see that Sa and 9 satisfy the following evolution equations [9]: 

{Sa)t = Va-9aU, (18) 

9t = ^ + ^. (19) 

Given Sa and 9, the curve T can be reconstructed up to a translation by integrating (I16p . 
However, we also need a point on the boundary to provide the constant of integration. 

Using the Sa — 9 formulation, we can reformulate the immersed boundary problem as 
follows: 

= -Vp + ^Au + L(X)(F(s„,0)), (20) 

V-u = (21) 
U = L*(X)(u(x))-n, (22) 
V = L*(X)(u(x)).T, (23) 

{Sa)t = Va-9aU, (24) 

9, = ^ + ^, (25) 



where 



F{Sa, 0) = ^{Tr) = Sb {Sa,aT + (s„ - l)9an) . (26) 



4 Spatial Discretization 

We use the spectral method to discretize both the Stokes equations and the immersed boundary 
equations in space since we work on periodic domains. We first discuss the discretization of the 
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Stokes equations in a regular N x N Cartesian grid with a uniform meshsize h. Let xj = jh 
and Uj = jh. The discrete Fourier transform and inverse Fourier transform are defined as 
follows: 

1 

J^h,.mk,y) = -Y,<p{x,,y)e-'^-^ =(t>{Ky), -N /2 + I < k < N /2, (27) 

3=0 
N-1 

J^h,ymx,k) = Y.'t'i^^Vj)^"''' =$i^^k), -N/2 + l<k<N/2, (28) 

j=0 

N/2 

^h'Mi^j^y) = E 0(fe,y)e*'^^^ =0(x,-,y), 0<i<iV-l, (29) 

k=-N/2+l 
N/2 

^h,li$)i^,yj) = E kx,k)e'''y^ =4>{x,yj), 0<j<N-l. (30) 

k=^N/2+l 

Now we introduce the discrete differential operator using the discrete Fourier transform de- 
fined above. For a function cpi^x, y) defined in the fluid domain Q, we approximate its spatial 
derivatives as follows: 

iDh,.(t>) (x, y) = T^l (ik (k, y)) , (31) 

{Dh,y^) (x, y) = T^l (ik {H,y<f) {x, k)) . (32) 

Denote V/^ = (Dh,a) ^h,y)- The differential operators are discretized in terms of D/^: 

Vp ^ VhP, (33) 
V-u ^ Vft-u, (34) 
V^n Vh-VhU = Vlu. (35) 

Next, we describe the discretization of the immersed boundary. We employ a Lagrangian 
grid with grid space Aa. The number of grid points along the boundary is Nf,. For a function 
ip{a) defined on the interface F, we define the discrete Fourier transform and its inverse as 
follows: 

J^Aamk) = — J2 H»j)e-""'' =i^{k), aj=jAa, (36) 

-^A^W(«.) = E me'"^^ =Hc^j). (37) 

When the interface is a closed curve, we can approximate the derivative operator along the 
interface by the spectral derivative: 

{DAa^) (a) = J'Zl {ik iJ^AaCl^) (k)) . (38) 

When the solution is not periodic, we can also use a finite difference method to discretize the 
derivative, we refer to ^27j for more details. 
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Now we discuss the discretization of the spreading and interpolation operators. These two 
operators both involve the use of a discrete delta function. The discrete delta function we use 
is introduced by Peskin in [27]: 

and 

{I (3- 2|r| + y/l + i\r\ -4r^) , |r| < 1, 

I (5 - 2\r\ - v/-7 + 12|r| - 4r^) , 1 < |r| < 2, (40) 

0, \r\ > 2. 

Using the above discrete delta function, we can discretize the spreading and interpolation 
operator as follows 

L/,(X)(g(a))(x) = ^ g(a)5;,(x - X(a, t))Aa, (41) 

LUX){u{^)){a) = ^ n(x)5;,(x-X(a,t))/i2. (42) 

The summation above is over grid points in T in (j4ip and over grid points in 0, in (j42p . Operator 
Lfi and are still adjoint using the following discrete inner product: 

</,<?>r^= f{a)9{a)Aa, (43) 
<u,v>Q^= u{x)v{x)h'^. (44) 

Using the inner product defined above, we have: 

< u(x),L(X)((7(a)) >n^ 
= Y n(x)L(X)(5(a))/i2 

= Y uixjh^ Y gia)Shix-X{a,t))Aa 

= Y 9{a)Aa Y ^i(x)(5h(x - X(a,t))/i^ 

= <L^(X)(«(x)),5(a)>r, . (45) 

As we will see later, this discrete self-adjoint property is crucial in obtaining our unconditional 
stable semi-discrete scheme for the immersed boundary problem. 

5 Steady Stokes flow 
5.1 Formulation 

For simplicity, we study the steady Stokes flow first. The governing equations for the steady 
Stokes flow are given as follows: 

= -Vp + /iAu + L(X)(F(s„,0)), (46) 
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V-u = 0, (47) 

U = u(X(a,t),t) -n, (48) 

V = u(X(a,t),t) -r, (49) 

{Sa)t = Va-9aU, (50) 

= E^ + Y^. (51) 

See 

In this simple case, we can use a boundary integral method for the two dimension Stokes 
flow (page 60 of [28]) to solve equations (P6j) - (|17j) to get the velocity on the boundary: 

«(X(a,t)) = [ l-{lnr)Fi{a',t) + ^Fi{a',t) + '^F2{a',t)]da', (52) 

viX{a,t)) = -^J^(^-{lnr)F2{a',t) + ^F2ia',t) + ^Fiia',t)^da', (53) 
where r = |r| and 

r = (ri,r2) = X(a,t)-X(a',t), F = (Fi,F2), u = {u,v). (54) 



5.2 Small Scale Decomposition 

As we can see from (j52p -(|53]). the velocity field can be expressed as a singular integral with a 
kernel ln(r). However, the singular velocity integral is nonlinear and nonlocal. It is difficult 
to solve for the implicit solution if we treat the velocity integral fully implicitly. The main 
idea of the Small Scale Decomposition technique introduced in [9] is to decompose the singular 
velocity integral into the sum of a linear singular operator which is a convolution operator 
independent of time t and the configuration of the curve, and a remainder operator which is 
regular. Since the remaining operator, which is nonlinear and nonlocal, is regular, the sim- 
plified convolution integral operator captures accurately the high frequency spectral property 
of the original velocity integral. Thus, if we treat only the leading order convolution operator 
implicitly, but keep the regular remainder operator explicitly, we can effectively remove the 
stiffness of the original velocity field which comes mainly from the high frequency modes of 
the solution. In this subsection, we will show how to perform such Small Scale Decomposition 
for the Immersed Boundary method applied to the Stokes flow. 

Observe that in the integral representation of the velocity field, (|52p - (j53p . the only singular 
part of the kernel is ln(r). The other part of the kernel is smooth. Thus to the leading order 
contribution of the velocity field can be expressed as follows: 

u(X(a,t)) — f -(lnr)F(a,t)da', (55) 

47r^ Jr 



1 



Sb 



(lnr)F(a',t) •T{a)da 
J — (Inr) {sa,a'T{a') + (sq — l)6a'n{a)) • T{a)da' . 



(56) 



(57) 
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Next, we perform a Taylor expansion for r, T(a') • T(a) and n(a') • T(a) as a function of a' 
around a. By keeping only the leading order term, we have 

r ~ SQ,(a)|a — a'l, (58) 
T(a') •T(a) ~ 1, (59) 
n(a') • r(a) ~ 0. (60) 

Substituting the above Taylor expansions to ([561) , we get 

47r/i 

Integrating by part, we obtain 



Y ~ / \i\.{Sa[0L)\a — o!\) Sa^a'd-ol . (61) 



where W is the Hilbert transform 

/(a' 



1 f(a') 

^[/](«) = - / ^^,dol. (63) 



-oo 

Using the same method, we can get the leading order contributions of U and XJ^ as follows: 

Sb 



u 



J -ln\a- a'\{sa' -l)Oa'da' (64) 

- -^/^^^^^f^cia' = -^W[(..-l)M. (65) 

Note that if / is a smooth function, then the commutator [7i,f]u = 7i{fu) — f'H{u) is a 
smoothing operator for u. Thus we can factor a smooth function from the Hilbert transform 
without changing its leading order spectral property. Suppose that Sa is smooth, then we 
obtain to the leading order that 

~ -j^{sa-i)n[e^]. (66) 

4// 

Applying the same analysis to the Eqs (|115p - (|116p gives 

(a A = -^WIDa^s^I + (da„1/ - DA^eU + ^W[DaoS«i) , (67) 
= 4 /! _ 1) „[Oa„«1 + + ^ 3 (l - ±) „|z,A„«l) . (68) 

4// V Sq,/ V Sa Sq, 4^ \ SaJ J 

Note that the leading order operator is linear. This suggests a natural semi-implicit discretiza- 
tion of the immersed boundary problem. 

Since we are dealing with a closed immersed boundary, it is natural to work in the Fourier 
space. Furthermore, the Hilbert operator has a very simple kernel under the Fourier trans- 
formation. Notice that 6 is not a periodic function of a. Its value increases 2tt every time a 
increases Lf,. Nevertheless, if we let 

27r 

0{a, t) = —a + (j){a, t), a G [0, Lb], (69) 
Lb 
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then (j) is periodic. It is more convenient to work with <p than 9. Substituting (j69p into (j68p 
and taking the Fourier transform on both sides of (j67p . (j68p . we obtain 



Sa,t 



5fe 
'4m 



-i\k\4> + 



(70) 
(71) 



where 7 = max 1 ) . We have also used the fact that Tik = —i sgn(/c) with sgn(/c) being 

" V So, J 

the signature function. The first term on the right hand side captures the leading order high 
frequency contribution of the terms from the right hand side. An important property of this 
leading order term is that it is linear in S(y and and has constant coefficient in space. This 
provides a straightforward application of the implicit time discretization. 

Since our small scale decomposition is exact near the equilibrium, we can use this result to 
get the stability constraint of the explicit scheme by using a frozen coefficient analysis. The 
stability constraint is given by 



At < C 



jjL h 
Sbl' 



(72) 



As we can see, the time step needs to be very small if Sfy is large and /i is small. For example, 
if Sb = 100, and ji = 10~^, then the stability would require that At < C10~'^/i. 



5.3 Semi-implicit schemes 

Based on the small scale decomposition presented in the previous subsection, we propose two 
types of semi-implicit schemes in this section. The ffist implicit time discretization uses the 
backward Euler method to discretize the leading order term while keeping the lower order term 
explicit. This gives rise to the following semi-implicit scheme: 



sn+l 



At 



+ 



in+1 



At 



4/i 



+ 



\ „n+l ' n+1 



(73) 

(74) 



We call the above discretization the semi-implicit method. Near equilibrium, the stability 
constraint of this numerical method is At < C{Sb,fx), independent of the meshsize h. Since 
the small scale decomposition only captures the leading order contribution from the high 
frequency components, this method can not eliminate the effect of Sb and ^ completely. The 
coefficients Sb and fx can still affect the time stability through the low frequency components 
of the solution, which comes from the second term of the right hand side. In order to obtain a 
semi-implicit discretization with better stability property, we can incorporate the low frequency 
contribution from the second term in our implicit discretization. This scheme can be found in 
the appendix A. We call it the semi-implicit scheme of second kind. 

The accuracy of the semi-implicit schemes presented above is just first order. In order to 
get a high order time discretization, we can use the integral factor method. The integral factor 
method factors out the leading order linear term prior to time discretization. They usually 
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provide stable and high order time integration methods for stiff problems. To use the integral 
factor method, we rewrite ([70]) . ([7T]) as 

^(e''*^^) = exp(^^|fc|t)p(s,,,^), (75) 

^(e€*<A) = exp(^^7|A:|t)g(l„,0), (76) 

where 

p(s„,« = ]^ {y„ - e„u) + ^\k\s„, (78) 

Q{s^A) = + (79) 

Now it is straightforward to discretize this system to high order. In particular, we can apply 
the classical fourth order Runge-Kutta method to discretize the above system to obtain a 
fourth order semi-implicit scheme. 

We remark that although the fourth order semi-implicit scheme based on the integral factor 
approach is much more accurate than the first order semi-implicit discretization, the stability of 
the fourth order method is weaker than the first semi-implicit scheme based on the backward 
Euler discretization. The fact that the higher order discretization gives a weaker stability 
property is a phenomenon which has been observed for almost all time integration methods. 
It is not a restriction of our semi-implicit schemes for the immersed boundary problem. 

The semi-implicit schemes we describe above only update the 9 and Sa variables. We also 
need to reconstruct the boundary at t""*"^ from ^"+^ and s^^. For this purpose, we need 
to update a reference point of the boundary. This will be done by using an explicit time 
integration method. The simplest one is the forward Euler method: 

a;"+i(0) = x"(0) + At(y"cos(r(0))-[/"sin(r(0))), (80) 
y"+i(0) = y"(0) + Ai(F"sin(r(0)) + ?7"cos(r(0))), (81) 

where U and V are evaluated at the reference point. A higher order integration method can 
be also used. In the explicit update of the reference point, we can use the values of U and V 
obtained using the semi-implicit discretization from the previous time steps to extrapolate the 
values of U and V in the intermediate time steps in our explicit update of the reference point. 
Once we have updated the reference point, we can obtain the configuration of the boundary 
(x,y) from {sa,0) by integrating (fT6]) 

a;"+i(a) = x"+i(0)+ r s'^+\a')cos{e''+\a'))da', (82) 

Jo 

y«+i(a) = y"+i(0) + rsS+H«')sm(r+i(Q'))rfa'. (83) 

Jo 

We can use more than one reference point, then average them to get the last configuration. 
This can improve the stability constraint significantly. Actually, in our computation, we use 
two reference points X{0),X{Ni)/2), then take the average to determine the position of the 
interface at next time step. Since we update only two reference points, the extra cost in 
updating the reference point is small compared with the overall computational cost. 
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6 Unsteady Stokes flow 



6.1 Formulation 

In this section, we will extend the semi-implicit discretization developed for the steady Stokes 
flow to the unsteady Stokes flow. The governing equations of the immersed boundary method 
for the unsteady Stokes flow are as follows: 

P-^ = -Vp + ^Au + L(X)(F(s„,0)), (84) 

V-u = (85) 
U = u(X(a,t),t) -n, (86) 
V = u(X(a,t),t) -T, (87) 

Sat = Va-e^U, (88) 

e, = ^ + ]^. (89) 

It is much more difficult to solve the fluid velocity u analytically from (|84 p - (j85p . As for the 
steady Stokes flow, we will first derive an unconditionally stable time discretization which will 
be given in next section and then apply the Small Scale Decomposition to the unconditionally 
stable time discretization to obtain our efficient semi-implicit schemes. 

6.2 An unconditionally stable semi-implicit discretization 

In this section, we will describe our unconditionally stable semi-implicit discretization of the 
Immersed Boundary method for the incompressible unsteady Stokes equations and prove its 
unconditional stability in the sense of total energy is non-increasing. 

The unconditionally stable semi- implicit discretization is consisted of two steps. In the first 
step, we update Sq,,u from to f"^^, then we get 0""''^ in the second step. 

Step 1: Update of u"+i and s^^^ . 

^ „n+^_ „n ^ _v^^"+i+^v2u"+i+L;,,„(F(s^+i,r;T",n")), (90) 

Vip"+^ = V,-L,,„(F(sri,r;T",n")), (91) 
= L^,„(u"+i).T" (92) 
= L^,„(u"+^).n", (93) 

= D^^V"^^ - D^^e^U^+\ (94) 



n+1 _ 



At 

where r" = (cos(r ), sin(0")), n" = (- sin(e"), cos(0")), Lh,n = L/,(X"),L;;„ = L^(X"), 
V/i and D^a are discrete derivative operators for the Eulerian grid and the Lagrangian grid 
respectively, and 

F(s2+\e";r",n") = 5^ (DAaS^'^" + (s^' " l)^Aa^"n") . (95) 
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Step 2: Update of 0"+^. After we have obtained u"+^, p"-~^^ and s^"*"^, we update 9 at 
t"""*"^ using the following semi-implicit scheme: 



V;,p"+i +/iV2u"+i (F(s-+i,r+i;T-,n") , (96) 



At 

V|p"+^ = V,-L,,„(F(.r\^"+';^",n")), (97) 



where 



= L^,„(u"+i).r" (98) 
17"+' = L;;,„(u"+i).n", (99) 

= ^ I)A.I7"^'+I)A.rF"+' . (100) 



F(sS+\r+^T",n") = 5fe (Da„sS+V" + (s^+i - l)Z5Aa^"+'n") (101) 



It is important to note that the above discretization is not fully implicit. In fact, both the 
spreading and interpolation operators are evaluated at the interface from the previous time 
step. Moreover, when solve the sJ^J"^^ and u"+^, in (j90p - (j94p . we use O"- instead of 0^~^^ to 
evaluate the force density. This makes our semi-implicit discretization linear with respect to 
the implicit solution variables, u"+^, 0"+^, and s^~^^. The above semi-implicit discretization 
essentially decouples the stiffness induced by the elastic force from the fluid equations. This 
enables us to remove the stiffness of the Immersed Boundary method effectively by applying 
the Small Scale Decomposition and arclength/tangent angle formulation as was done in [9]. 

In the following, we will prove that this semi-implicit discretization is unconditionally stable 
in the energy norm. 

By using a discrete summation by parts, we can show that 

< f,DAag >ri,= - < DAaf,g >ri,, < n,Vhg >n^= - <Vh ■ ii,g >Qi^ . (102) 

First, we define the total energy of the physical system. The total energy includes the 
kinetic energy K and the potential energy P, which are defined below: 

1 ^ 

1 S" 
P = -55<s,-l,s,-l>r,= yE(«"J-l)'^«- (104) 

i=i 

The total energy is then defined as 

E = K + P. (105) 

Below we will prove the unconditional stability of our semi-implicit discretization. To 
simplify the presentation, we still denote the discrete spectral derivative D/^aQ of a function g 
as Qa- 
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Taking the discrete inner product defined by dM]) of with u"+^ + u" and using (jl02p . 
we obtain 

2(i^«+i -K"") = p < u"+^ + u'^, u"+^ - u" >n^ 

= p< -u"+i + u", u"+i - u" >f,, +2p < u"+i, u"+^ - u" >n, 
= -/9<u"+i-u",u"+i-u">n, 

+2At{< u"+\ -V,p + /iV^u^+i + L,,„ (F(.r\ r)) >n J 
= -p< u"+i - u",u-+i - >n, -2At < u^+\VhP >n, 

+2At < u"+i,^v2u"+i >n, +2At < u"+\L,,„ (F(sr\r)) 
= -p< u"+i - u",u"+i - u" >n, -2At < Vh ■ u"+\p >n, 

< V;,u"+\V;,u"+l 

+2At < L^^ (u"+i) , FK+\en >r, • (106) 



The second term on the right hand side of (jl06p is zero because the discrete velocity field is 
divergence free, i.e. V/i • u"^^ = 0. The fourth term can be rewritten as 

<Ll^, (u"+i),F(.r\n >r, 
= < + C/"+in", Sb (sS+V" + (s^+i - 1)0>") >r, 

= 5,(<F"+\.S+i>r, + <C/"+\(.r^-l)C>r,). (107) 
Combining (jl06p and (jl07p . we can get 

2(i^"+i - i^") = -p< u"+i - u", u'^+i - -2pAt < V;,u"+\ V;,u"+^ 

+2S,At{< V^+\sl+^ >r, + < C/"+\(5r^ - 1)C >rj. (108) 

Similarly, by taking the discrete inner product defined by (03]) of ([MI) with 5^^+^ + - 2 and 
using p02p . we get 

2(pn+l_pn) ^ Sb<sl+^ + sl-2,sl^^-sl>T, 

= S,< -sl+' + sl, sl+' - sl >r, +2S, < - 1, - >r, 
+25fcAt < - 1, V^+^ - >r, 

+25feAt(- < s^%\ >r, - < - 1, >r J. (109) 

Adding (fTOSll to (fT09]) . we have 



1 

< 0. (110) 
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This proves that our semi-imphcit discretization is unconditionally stable in the sense that the 
total energy is non-increasing. 

Remark 1. In our proof presented above, we have used two important properties of our 
semi-implicit discretization. The first property is that the discrete spreading and interpolation 
operators are adjoint. The second property is that the velocity field satisfies the discrete diver- 
gence free condition. It is clear from the above proof that as long as these two properties are 
satisfied by our spatial discretization, the corresponding semi-implicit discretization introduced 
in the previous subsection is unconditionally stable. 

Remark 2. We remark that the proof above is similar in spirit to that of a semi-linear dis- 
cretization obtained by Newrenn et al in [25]. There is some minor difference between the 
unconditionally stable semi-implicit discretization obtained by Newren et al and our uncondi- 
tionally stable semi-implicit discretization. In the problem considered by Newren et al., the 
force is a linear function of the interface. On the other hand, in the problem we consider, the 
force is a nonlinear function of the interface (the rest length of the boundary is not zero). By 
using the Sa — formulation, the force is a linear function of Sa- By treating 9 explicitly, we 
obtain a semi-implicit discretization that is linear with respect to Sq. Due to the decoupling 
between Sa and 9, we need to solve two Nf, x A';, linear systems instead of one 2Niy x 2Nij linear 
system in the semi-implicit discretization obtained by Newren et al. 

Remark 3. For the steady Stokes flow, we can also prove the following semi-implicit dis- 
cretization is unconditionally stable: 



Step 1: 



= -V,p'^+^ + /.V2u"+^ + L,,„(F(sr\r)), (111) 

V|p'^+^ = V,-L,,„(F(.ri,r;r",n"))), (112) 

= L;;,„(u"+1(x)).t«, (113) 

= Lt„(u"+i(x)).n", (114) 



„n+l _ r. 
°a ^ 

At 

0n+l _ QT 
At 



DAaV''+^ - DAa9''U''+\ (115) 
i (DAaU''+' + DA,r+V"+M . (116) 



Step 2: 












n+1 


= V;,-L;,,„(f(s^+1,i 


V 


n+1 


= L^,Ju"+i).r" 


u 


n+1 


= L^,Ju"+i).n", 



(117) 
(118) 
(119) 
(120) 

1 / -n+1 



At 



^ {DAaU''^' + DA^r y . (121) 
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In this case the total energy is just the potential energy 

E=^<s^-hsa-l >r,= y - l)'Aa. (122) 

i=i 

As in the case of the unsteady Stokes flow, as long as the velocity field satisfies the discrete 
divergence free condition and the discrete spreading and interpolation operators are adjoint, 
we can prove that above semi-implicit discretization is unconditionally stable in the sense of 
total energy is non-increasing. 



6.3 Small Scale Decomposition 

In order to apply the Small Scale Decomposition to our unconditionally stable time discretiza- 
tion, we would like to solve for the velocity field at time t"^^ from the space-continuous version 
of (|90p and (j9ip using an integral representation: 



u"+^(x) = (1-^V^j (^u" + — (1- V(V^)-^V-)Ln(F(s^+\r)) 

i-^v^)"' fu«+^L„(F(.r\n) 

p J \ p 

^ f 1 _ /^v^V' (V^)-! (VV • L„(F(.ri,r)) 



P \ P 

To solve for the velocity field at t""*"^, we need to use the following free space fundamental 
solutions in two space dimensions which are defined as follows: 



P 



V^^ Ei = 6{x-x'), (123) 



V2 (^1 _ j E2 = <5(x - x'). (124) 

These two fundamental solutions can be expressed in terms of the modified Bessel function of 
the second kind [1]: 

El = ^i^o(Alx-x'l), (125) 
E2 = ^(i^o(A|x-x'|+ln(|x-x'|)), (126) 

where = ^ and Kq is a modified Bessel function of the second kind. By integrating by 
part, we can further express the velocity u""*'^ as 

u"+i(x) = / A2i^o(A|x-x'|)u"(x')dx' 
27r Jn 

+7^— / A2i^o(A|x-X"(a')|)F(.r\n^«' 

I VxVx(i^o(A|x-x'|) + ln(|x-x'|)).L"(F(.ri,r))'ix 
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1 At 

where G is defined as follows: 



^ I A2/^o(A|x-x'|)u"(x')<ix' 

I A2i^o(A|x-X"(a')l)F(sr\n^«' 

f G(x-X"(a')) •F(s^+\r)da', (127) 
Jn 



G.,{r) = ^-^ + ^A^(i^o(A|r|) + K2(A|r|))p 

-Ai^i(A|r|)(|i-p), (128) 

and Kq, Ki, K2 are all modified Bessel functions of the second kind [1]. 

In this subsection, we will perform a Small Scale Decomposition to the velocity field based 
on the integral representation (jl27p . Recall that in our semi-implicit discretization, the velocity 
field at t""*"^ is evaluated on the boundary X" at t". Thus we should perform our Small Scale 
Decomposition for u"+^(X"). To this end, we first write down the integral expression of 
u"+i(X") as follows: 



1 At 
'^2^~p 
1 At 



u"+i(X"(a)) = — / A2i^o(A|X"(a) -x'|)u"(x')(ix' 
2vr jQ 

■ I A2i^o(A|X"(a) - X^{a')\)¥{sl+^,e'')da' 

■ / G(X"(a) - X'^(a')) • F(s!:+\ ^")da'. (129) 
27r p Jn 

To perform the Small Scale Decomposition to the above velocity integral, we would like to 
decompose the singular velocity kernel as the sum of a linear singular operator of convolution 
type and a remainder operator which is regular. Using the Taylor expansion for a' around a, 
we get the following decomposition: 

T^"+^(a) =u"+^(X"(a)) •T"(a) 
SbAt 



27rp 
SbAt 



f X'Ko{Xs"Ja-a'\)sl\]da' 

. / (lx\Ko{Xsl\a-a'\)+K2{Xsl\a-a'\)) 
I-irp Jt" \2 (sjj) (a - a') 



s 



n+l 



]da', (130) 



where inside KQ{Xs'^\a — a'\) is evaluated at a. Notice that [T] 

^ ri^o(AsSla-a'l) ) = jA2(Ko(As2|a-a'l) + ^2(As^|a-a'|)). (131) 



Integrating the right hand side of (jl30p by parts twice, we get 

^"+'(«) ~ i^SbAt [ X^KoiXs^'Ja-a'DsllUa' - 
Inp JT" ' 

/ {K,{Xsl\a -a'\)- ln(a - a')) sl'^^,^,^,da' . (132) 
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Similarly, we can obtain the leading order contribution of U "^^ as follows: 

ShAt 



[ {KoiXs^Ja - a'\) - ln{a - a')) - 1)61+') , , da'. (133) 
Using this decomposition, we obtain the following scheme: 

T{s^+') + (DAaV*'^+' - i^Aae"f/*'"+' - T{s^)) , (134) 

-V,j9'^+i + ^V2u"+i + L,,„(F(.ri,n), (135) 

V^p"+i = V,,-L;,,„(F(.r',^")), (136) 

= Ll^{u^+')-T^, (137) 

= Ll^{u^+')-n^, (138) 

^r!.+l_5|n ^ 5(^"'+^) ^ f 1 ^^^+1 ^ ^ an^rn+l\ aian\\ ^^39) 



„n+l _ „n 



where 



-^^/ (i^o(AsSI«-«'l)-ln(a-a'))CWrf«'l > 



\2np{s[ 

5(^"+') = f ^ ^'fL / (i^o(A.S|a-a'|)-ln(a-a'))((.r^-l)^a') , , 
and u*'"+^ is the velocity at which is calculated explicitly 

P = -V,/'"+i + /.V2u-"+i + L,,„(F(.^,n) (140) 

V2p*'"+i = V,-LUnSa,01) (141) 
^*,n+i ^ L^,^(u*'"+i) • r" (142) 
f^*,n+i ^ l;;^^(u*'"+^) -n". (143) 

The derivation of the above semi- implicit scheme is given in Appendix B. 

However, the expressions of T and S are still too complicated and need to be further 
simplified. The leading order linear operator, which contains K(){Xs^{a)\a — a'\), is not a 
convolution operator. Thus, it does not have a simple kernel under the Fourier transform as 
the Hilbcrt operator in the case of the steady Stokes flow. To further simplify the kernel, 
we approximate s]l(a) by min s"(a). With this approximation, the corresponding leading 
order operator is a convolution operator and can be diagonalized under the Fourier transform. 
Denote /3 = AminsQ(Q;). In Appendix C, we will show that 
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Using (jl44p and replacing s^{a) by mins^(a), we can simplify the leading order term T{s'2^~^^) 
and S{6^^^) under the Fourier transform: 



n+l^ 



ShAt 



s{e 



n+l\ 



5fcAtmax(sS+^ - l) ( 



mm s. 



2p (min s"^^ ' 



V 



(^AminsJJ^ + /c^ 



(145) 



(146) 



1 



When /X » 1, we have A = ^ 1. By Taylor expanding (11450 and (I146P with respect to 

A and keeping only the first order term, we obtain the leading order term as follows: 



5(r+i) 



max 



1 



(147) 
(148) 



which is the same as the steady Stoke flow. This is also consistent with one's physical intuition. 
When the viscosity is very large, the flow changes very slowly. The inertial term can be 
neglected. 



When <C 1, then A 



1 



f(.ri) 



^ 1, the asymptotic expansion is 



2gn+l 



2 ( min si 



-k's 



S6Atmax(s^+^ - l) 



3fln+l 



2p I min 



kW 



(149) 
(150) 



From the asymptotic expansion above, we can see that our small scale decomposition is also 
consistent with the linearized stability analysis which Stockie and Wetton got in [3^. Using 
the leading order term above, we can get the leading order term of the eigenvalue same with 
the result in [30j . 

We can also obtain the corresponding stability constraint for the explicit scheme near the 
equilibrium: 



At < C{Sb,p)h^, 



(151) 



where 1 < /3 < 3/2. The value of (3 depends on /i. li p 1, then we have (3 ~ 3/2. On the 
other hand, if // S> 1, we have /3 « 1. 
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6.4 The numerical scheme 



Based on the small scale decomposition we developed in the last subsection, we can now 
describe our semi- implicit numerical scheme. Combining the time discretization (j90p - (jlOOp 
with the decomposition (|130p - (jl33p and using the approximation (jl45p - (jl46p . we obtain the 
following semi-implicit numerical scheme: 



Step 1: Update of u'^+i and s^-^^. 



At 



u 



n+l 



U 



At 



where 



n+l\ 



S.At ( (Xmins'^Ye + k^ \ 



2p (mm Mxmms^y + 



(152) 

(153) 
(154) 



(155) 



(156) 



and u*'""^^ is the intermediate velocity at t'^'^^ which is calculated by solving the unsteady 
Stokes equations implicitly while evaluating the elastic force explicitly: 



P- 



72, *,n+l 



At 



V,-L,,„(F(sS,r)), 

y*,n+i ^ l;;^^(u*'"+1) • r", 

TT*,n+l T* /',,*, ™ + 



n . 



(157) 

(158) 
(159) 
(160) 



Step 2: Update of 6'"'+^. Once we have updated u, p, and at t"+^, we update Q 



n+l 



3n+l 



using the following semi-implicit scheme: 



where 



At 



n+l^ 



5(0 



+ 



r* 



S'ftAt max (s^ 



2/0 ( 



mms. 



5(r) 



V 



(a nnn s"^ + A;^ ^ 



3n+l 



(161) 

(162) 
(163) 

(164) 



This is our semi-implicit scheme for the unsteady Stokes flow. The spectral discretization 
in space has the advantage of being high order accurate and the leading order operator has 
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a simple kernel under the Fourier transform. As it is, the time discretization is only first 
order. Based on the first order semi-implicit scheme that we develop in this subsection, we 
will develop a second order semi-implicit scheme in the next subsection. 

A near equilibrium stability analysis shows that the stability constraint of this semi-implicit 
scheme is of the form At < C{Sb,fi), which is independent of the wave number, but still 
dependent on Sh and n. This is due to the fact that the Small Scale Decomposition does 
not capture the low frequency components of the solution accurately. The low frequency 
components of the solution can affect the stability of the time discretization in two ways. 
The first one is through the small scale decomposition, which only captures the leading order 
contribution of the solution at high wave numbers. The second one comes from the second 
term of the right hand side of the dynamic equations for Sa and 6. As in the case of the steady 
Stokes flow, we can include the leading order contribution from the second term in our leading 
order term and treat them implicitly. This treatment would significantly improve the stability 
property especially when the elastic coefficient is large or the viscosity is small. This improved 
stability is at the expense of solving a linear system for the implicit solution at each time step. 
We call this semi-implicit discretization as the semi-implicit method of the second kind. More 
discussions on the semi-implicit method of the second kind can be found in Appendix A. 

Remark 4. The leading order term we derive above is calculated analytically using the space- 
continuous formulation with an unsmoothed Dirac delta function. As Stockie and Wetton 
pointed out in [32j , this analysis over-predicts the stiffness of the Immersed Boundary method 
in a practical computation. If we use the leading order approximation directly, the semi-implicit 
scheme with the leading order terms derived above tends to over-dissipate the solution. To 
alleviate this effect in the practical implementation, we rescale the leading order term by a 
coefficient which is calculated at the first time step in the following way: 

maxo V^'* 
maxa T (sO ) ' 

maxa 
maxc Sjj (9^) ' 

where Su{6^) is the leading order term of U ^ , which can be computed from S{9^) via the 
Fourier transform. The leading order term we use in a practical computation is actually 
C7vr(s^+i) and CuS{9''+^). 

6.5 A second order semi-implicit scheme 

Based on the first order semi-implicit scheme we have developed in the previous subsection, 
we will derive the corresponding second order semi-implicit scheme in this subsection. 

First, we need to use a robust implicit second order temporal discretization. To simplify 
the presentation, we will only describe the semi-discrete algorithm. The space discretization 
is done in the same way as before. The second order temporal discretization we use consists 
of two steps. In the first step, we take a fractional time step from to t"^2. It is same with 
the first order semi- implicit discretization (f90]) - (jl00p , except the timestep is 

In the second step, we integrate the unsteady Stokes equations from to t^~^^ based on 
the midpoint and the trapezoidal rules: 



Cv = 
Cu = 
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step 1: Update of u"+^ and sg+^ 



p- 



= -Vp + /iV2u + L„+i(F(s-„,r+5;r"+in»+i)), (165) 
V^p = V-L„^i (F(s„,r+^;T"+in-+5)), (166) 



V = L;^i(u)-T"-^2, (167) 



2 



= L;+i(u)-n"+5, (168) 

At/2 (169) 

where t"+5 = (cos(r+5), sin(r+i)), n"+5 = (- sm(^"+5), cos(^"+5)), L„^i = L(X" + 
i),L;^i =L*(X- + i)and 

F(Sa,r+5;r"+5,n"+^) = 55(Z?Aa5aT"+^+(^a-l)i^Aar+5n"+i) (170) 

Step 2: Update of 0"'+^. After we have obtained u'^"^^ and s^'^'^, we update 9 at t""*"^ using 
the following semi-implicit scheme: 



n+l _ „n /fjn+1 1 „n+l\ _ 

= -Vp'^+i i +L_i (F(s„,^";T"+2,n"+2)) , (171) 



At 

V^p-+^ = V-L„+i (F(sa,^;T-+^,n-+^)), (172) 



Qn+1 _ Qn I 



At 

where 
and 



-iUc^ + ea'V]. (175) 



F(sa,^;T"+in"+5) = 56 (s„,aT"+^ + (sa - 1)^,11-+^) (176) 



u = ^ , ga= " ^ ^ = ^ ■ (177) 

Here, L^_^_i and are the spreading and the interpolation operators evaluated at X"''''2. 

Using the same method of analysis, we can prove that the above second order semi-implicit 
discretization is unconditionally stable in the sense that the total energy is non-increasing. 

The first step in our second order method is identical to the first order method except 
that the time step is instead of At. Thus we can use the first order semi-implicit scheme 
introduced in last subsection to compute it directly. 
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In the second step, we can also apply the Small Scale Decomposition with some modi- 
fications. After applying the Small Scale Decomposition to the second step of the two-step 
method, the second step of the semi-implicit scheme has the form 



Step 1: Update u''+^ and 5^+^ 



n+l _ n 
■'a rp 



At 



= - Vp + nV^u + L , 1 F s, 



" At 

= V.L„^. (F(s-e.,r+^)), 

The leading order terms, T is given by 

/ 

SbAt 



(178) 

(179) 

(180) 



4p mm Sa 



\ ■ ""'"2 

A mm Sa 



V 



A min ^ ) -|- A;^ 
a J 



\kf 



(181) 



(182) 



where A 



2p_ 



, and u*'"+^ is the intermediate velocity at t""*"^ which is obtained by solving 



the unsteady Stokes equations implicitly but with the forcing evaluated explicitly: 



At 



V* 

u* 



-vr (f , 
v-L„+i (f 



= L* i(u*)-n"+i 



and 



u 



*,n+l 



U 



(183) 
(184) 
(185) 
(186) 

(187) 



Step 2: After the u"+^ and s^^^ are calculated in step 1, we update 6 at t"'^^ using the 
following semi-implicit scheme: 

-] + 



^n+l _ Qn 

At 



= S 



-\+^(Ua + eaV-s{e'^-'^ 

' Srv 



(188) 



where 



V 



T* 



U = L* 1 

n+2 



(189) 
(190) 
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and the leading order term is 



) 



/ 



Si,At max [sa ^ — 1 



m = - 




\k\ 



3 





9 



(191) 



4p I mm Sa 



\ 



A min s, 



a 



I 



This completely defines our second order semi-implicit scheme. 

7 Numerical results 

In this section, we will perform a number of numerical experiments to test the stability of our 
semi-implicit schemes for both the steady and unsteady Stokes equations. We also compare 
the performance of our semi-implicit schemes with the explicit scheme and the fully implicit 
scheme. Our numerical results indicate convincingly that our semi-implicit schemes has a much 
better stability property that that of the explicit scheme. Moreover, the computational cost of 
our semi-implicit schemes is comparable to that of an explicit scheme. Thus our semi-implicit 
schemes offer significant computational saving over the explicit scheme, especially when the 
number of grid points is large. 

7.1 Model problem 

The test problem we use is one typically seen in the literature, in which the immersed boundary 
is a closed loop initially in the shape of an ellipse. We choose an ellipse initially aligned in the 
coordinate directions with horizontal semi-axis a = 0.32 and vertical semi-axis b = 0.24. The 
boundary can be parameterized as follows: 



The fluid is initially at rest in a periodic domain O = [0, 1] x [0, 1]. We use a periodic boundary 
condition for the fluid flow. For the initial condition defined above, the rest state of the 
boundary is a circle with radius r = 0.2. For the unsteady Stokes flow, the immersed boundary 
with the above initial condition evolves as damped oscillations around a circular equilibrium 
state. The area is conserved during the time evolution since the flow is incompressible. For 
the steady Stokes flow, the boundary converges to the circular state without oscillations. 

We use a uniform N x N grid to discretize the fluid domain, J7. We choose Nj, = 2N 
number of grid points to discretize the immersed boundary so that there are approximately 
2 immersed boundary points per mesh width. We use the spectral method to discretize the 
spatial derivatives both in the fluid domain and along the immersed boundary. The leading 
order singular integral is also discretized by the spectral method. 




(192) 
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7.2 Steady Stokes flow 

First, in order to reduce the number of parameters in our test problem, we write the equations 
in terms of the following dimensionless variables to get the nondimensional model [33j, 

., t , X , uto , pto „, fLto 

i = — , X = -, u = — , p= — , f= , 

iQ L L ^ fi 

where L is the size of computational domain, to is characteristic time. 
Using these new variables, we have 

= -Vp' + Au' + f'(x',t')> (193) 
= V-u'. (194) 

For the equations of the elastic boundary, the dimensionless variables are 

= J^^ «a = ^ = ^' " = 5^' ~si^' T = T". n = n. 

Then the equations describe the interaction of the boundary and the fluid become 

U' = u' (K'{a',t'),t') -n', (195) 
V' = u' {X'{a',t'),t') -t', (196) 
s'a,t' = Vi>-e',,U', (197) 

^' - ^K'+V'0'^,), (198) 



't' 



where 



f'(x',t') = /■ " F'{a',t')6{x' -X'{a',t'))da', (199) 

fiL Jo 

u'(X'(a',0,0 = f u'{^',t')6{ji.' -X.'{a',t'))dji.'. (200) 
Jn 

There are two nondimensional parameters in this problem: 

fj,L ' L 

If we let to = thsn the only parameter in this dimensionless model is Li,/L which is fixed 

in our test problem. So we can always fix S";, = /i = 1 in our numerical study. 

The stability analysis in the steady Stokes flow suggests us to use the total energy as a 
criterion to test the stability of different numerical methods. For the steady Stokes equations, 
the total energy is equal to the potential energy. In Fig[Tl we show that the energy for four 
different numerical methods: the explicit scheme, the semi-implicit scheme of first kind, the 
4th order semi-implicit scheme using the integral factor method and the unconditionally stable 
semi-implicit scheme. In this figure and the subsequent figures, we use the the legend "semi- 
implicit" to denote the semi-implicit scheme of first kind, and the legend "integral factor" to 



25 




Figure 1: Energy of the system for four different schemes. = 128, 5^ = 1, = 1. Left one: 
At = 0.1; Right: At = 1. Here the legend "stable semi-implicit" stands for the unconditionally 
stable semi-implicit scheme, "semi-implicit" for the semi-implicit scheme of first kind, and 
"integral factor" for the semi-implicit scheme based on the integral factor method. 



denote the semi-implicit scheme based on the integral factor method. We use two different 
time steps, 0.1 and 1, respectively. When At = 0.1, all the four methods are stable. They 
give almost identical results. When At = 1, the explicit scheme becomes unstable , but all 
the semi-implicit schemes are stable. In fact, all the semi-implicit schemes remain stable with 
much larger time steps. In Fig [2l we plot the energy of the system for semi-implicit schemes 
of first kind and the semi-implicit scheme based on the integral factor method with At = 10. 
Fig [3] shows the configuration obtained by the two semi- implicit schemes at the final time with 
At = 10. They both remain as a circle, but lose some area compared with the original state. 

Next, we compare the performance of our semi-implicit schemes with the explicit and fully 
implicit schemes. The fully implicit scheme we use here was originally proposed by Tu and 
Peskin in [33j. In order to make a fair comparison, we run the implicit schemes (semi-implicit 
and fully implicit) with a time step small enough to make sure that the computational results 
have a reasonable accuracy. We take At = 4 for the fully implicit and the semi-implicit 
schemes. With this time step, the area loss is less than 5%. For the explicit scheme, we 
take At = 1/4, 1/8, 1/16, 1/32 which corresponds to = 64, 128, 256, 512 respectively. These 
time steps are the largest possible to keep the stability of the explicit scheme. The time we 
compute is T = 20. The result is shown in Table [H From this comparison, we can see that the 
performance of our semi-implicit schemes is much better than the explicit, the fully implicit 
scheme, and the unconditionally stable semi-implicit scheme in all cases. As we can see, the 
larger the number of the spatial grid points is, the more computational saving we would get 
using our semi-implicit schemes. Even for a modest grid size, our semi-implicit schemes still 
give a significant computational saving compared with the explicit or the fully implicit scheme. 
It is interesting to note that although the computational cost of the unconditionally stable semi- 
implicit scheme (labeled as s,s,i) is faster than the fully implicit method, the computational 
cost of the unconditionally stable semi-implicit method is still more expensive than the explicit 
scheme. This makes the unconditionally stable semi-implicit scheme not very practical. 
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Figure 2: Energy for the semi-implicit scheme of first kind (labeled as "semi-implicit") and 
the semi-implicit scheme based on the integral factor method (labeled as "integral factor"). 
At = 10, N = 128, Sfc = 1, /X = 1. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 3: Dashed line: the initial boundary configuration; Solid line: the boundary configura- 
tion after 20 time steps with At = 10, A'' = 128, Si, = 1, /j, = 1. Left: semi-implicit ; Right: 
integral factor method. 
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N 


cxp 


s,i 1 


s,i 2 


s,s,i 


f,i 


64 


1 


0.4 


2 


7 


9 


128 


5 


0.7 


3 


30 


39 


256 


30 


1.3 


7 


139 


206 


512 


344 


4.2 


19 


611 


1200 



Table 1: Execution time for each computation in seconds. The legends are defined as follows: 
"exp" stands for the explicit scheme, "s,il" the semi-implicit method of first kind, "s,i2" 
the semi-implicit scheme of the second kind , "s,s,i" the unconditionally stable semi-implicit 
method, and "f,i" the fully implicit scheme . N is the number of grid points along each 
dimension. 

7.3 Unsteady Stokes flow 

We can also get the nondimensional model for unsteady stokes flow. Similar as the steady 
stokes case, we define the following dimensionless variables: 

, t , X , uto I pto , fLto 
iQ L L pi fj, 

where L is the size of computational domain, to is characteristic time. Using these new vari- 
ables, we have 

W = 0(-Vp' + Au' + f'(x',i')), (201) 

= V-u'. (202) 

For the equations of the elastic boundary, the dimensionless variables are 

X = — , s^ - —, d-d, O! - —, 1 - —, b - —, T -T, n - n. 

Then the equations describe the interaction of the boundary and the fluid become 

U' = u'{X.'{a',t'),t')-n', (203) 
V' = u'{X'{a',t'),t')-T', (204) 
Sa,t' = V!,,-e'^,U', (205) 

o't' = ^{K' + v'e'^,), (206) 



where 



f'(x',t') = F'{a',t')d{x' -X.'{a',t'))da', (207) 

fiL Jo 

u'{-X'{a',t'),t') = I u'(x',i')5(x'-X'(a',i'))(ix'. (208) 
Jo, 
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From the nondimensional analysis, we can see that there are three nondimensional parameters 
in this problem: 

Sbtp fitp Lb 

If we let to = ^-—■> then the parameters left in this dimensionless model is -^^-^ = — — and 
Sb pL^ pLSb 

Lb 

— . Lb and L are fixed in our test problem only depends on the initial condition. So — is 

L pLSb 
the only parameter in our test model. For this reason, we always fix the elastic coefficient Sb 

to 1, but vary p. 

In our computations, we use the following parameter values: 

p=l^Sb = l,p = 0.1, 0.01, 0.005. 

We vary the number of the spatial grid points along each dimension in the following fashion: 

A^ = 64, 128,256,512. 

The criterion that we use to check whether one scheme is stable or not is that the total 
energy of the system is non-increasing and the boundary configuration lies within the compu- 
tational domain. 

Next, we perform some numerical experiments to test the stability of our semi-implicit 
schemes for the unsteady Stokes flow. Fig. S] shows that the energy obtained by the explicit 
scheme and the semi-implicit scheme of the first kind. We take two different timesteps, 0.005 
and 0.05. With At = 0.005, the explicit and semi-implicit schemes are all stable, and they give 
nearly identical results. With At = 0.05, the explicit scheme becomes unstable, but the semi- 
implicit schemes remain stable. Even if we increase the timestep to At = 1, the semi-implicit 
methods are still stable, as we can see from Fig. [5j For the semi-implicit scheme of the first 
kind, we have used the Small Scale Decomposition and further simplification of the singular 
integral kernel. Therefore, the total energy in our semi-implicit scheme is not guaranteed to 
decrease monotonically in time. Nonetheless, we observe that the total energy still decreases 
in time as is the case for the unconditionally stable semi-implicit scheme. In Fig. [5l we also 
plot the boundary configuration at the final time step, which is an approximate circle. 

We remark that the semi-implicit scheme is not unconditionally stable, although its stability 
is much better than the explicit scheme. This is due to the fact that we have used the Small 
Scale Decomposition and further approximation of the leading order singular integral operator 
to simplify the computation of the implicit solution. As we mentioned before, the Small Scale 
Decomposition captures only the high frequency contribution to the stiffness, but it does not 
remove the stiffness of the system induced by the low frequency components of the solution. 
Thus there is still some mild timestep stability constraint for the semi-implicit scheme. The 
time step has a mild dependence on the elastic coefficient Sb and the viscous coefficient p. On 
the other hand, our numerical study shows that the time step for the semi- imp licit scheme is 
independent on the meshsize. 

We also compare the performance of our semi-implicit schemes with the explicit scheme. We 
do not compare the performance of our semi-implicit schemes with the fully implicit scheme 
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Figure 4: Total energy of the unsteady Stokes system for different schemes with two different 
timesteps. N = 128, Sb = l, n = 0.01. Left: At = 0.005; Right: At = 0.05. The legend "semi- 
implicit" stands for the solution obtained by the semi-implicit scheme of first kind, "stable 
semi-implicit" the unconditionally stable semi-implicit method. 




Figure 5: Total energy of the system for two semi-implicit schemes. Here "stable semi- 
implicit" stands for the unconditionally stable semi-implicit scheme, and "semi-implicit" the 
semi-implicit scheme of first kind, At = 1, N = 128, Sh = 1, n = 0.01. 
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Figure 6: Dashed line: the initial boundary configuration; Solid line: the boundary configura- 
tion after 20 time steps with At = 1, N = 128, Sb = 1, = 0.01. Left: the unconditionally 
stable semi-implicit scheme; Right: the semi- implicit scheme of the first kind. 

here because the fully implicit scheme is quite expensive and is not competitive with the 
explicit scheme. In order to keep the area loss is no more than 5%, we take At = j for all of 
the semi-implicit schemes . For the explicit scheme, we take At = 1/64,1/128,1/256,1/512 
which correspond to the spatial mesh sizes N = 64, 128, 256, 512 respectively, when /i = 0.05. 
When /i = 0.01 and 0.005, the time step is set to be At = 1/128,1/256,1/512,1/1024 and 
t = 1/256, 1/512, 1/1024, 1/2048. These time steps are the largest ones we can take to keep the 
stability. We compute the solution up to T = 2. The results are documented in Table [21 We 
can clearly see that the semi-implicit scheme of the first kind gives a significant improvement 
over the explicit scheme. The cost for the semi-implicit scheme of the second kind is higher 
than that for the semi-implicit scheme of the first kind. This is because we need to solve 
for a linear system at each time step for the semi-implicit scheme of the second kind. The 
cost increases as the number of the spatial grid points increases. The semi-implicit scheme 
of the second kind and the unconditionally stable semi-implicit scheme both need to solve a 
linear system at each time step. Their complexity are same, both are 0{N^). But for the 
unconditionally stable semi-implicit scheme, the scaling constant in front of A'^^ is much larger 
than the semi- imp licit scheme of the second kind. The reason is that the cost of computing 
the coefficient matrix of the linear system for the unconditionally stable semi-implicit scheme 
is much higher. As we can see from Table 2, the unconditionally stable semi- imp licit scheme 
(labeled as s,s,i) is still quite expensive compared with our semi-implicit schemes that use the 
Small Scale Decomposition. For fj, = 0.05, the unconditionally stable semi-implicit scheme is 
even more expensive compared with the explicit scheme. Although the unconditionally stable 
semi-implicit scheme is slightly faster than the explicit scheme for smaller /x, the semi-implicit 
scheme (labeled as s,i,l) which uses SSD to further simply the singular integral kernel gives a 
much more efficient algorithm. It gives a factor of 242 times speed-up over the explicit scheme 
in the case of /u = 0.005 and = 512. 
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N 


/i = 0.05 


/i = 0.01 


fi = 0.005 


exp 


s,i 1 


s,i 2 


s,s,i 


exp 


s,i 1 


s,i 2 


s,s,i 


exp 


s,i 1 


s,i 2 


s,s,i 


64 


1.8 


0.5 


4 


11 


3.3 


0.5 


4 


12 


6.6 


0.5 


4 


12 


128 


9 


1 


10 


48 


18 


0.9 


10 


47 


35 


0.9 


10 


48 


256 


58 


2.4 


25 


229 


116 


2.4 


25 


228 


236 


2.4 


25 


226 


512 


738 


12 


99 


980 


1461 


12 


98 


982 


2910 


12 


98 


977 



Table 2: Execution times for each computation in seconds. The legends are defined as follows: 
"exp" stands for the explicit scheme, "s,il" the semi-implicit scheme of the first kind, "s,i2" 
the semi-implicit scheme of the second kind, "s,s,i" the unconditionally stable semi-implicit 
scheme. 

7.4 Second order semi-implicit scheme for the unsteady Stokes flow 

In this subsection, we perform numerical experiments to test the convergence rate and the 
stability property of our second order semi-implicit scheme. To check the convergence rate in 
time, we set N = 256 and vary the time step in powers of 2 from ^ to j^. When = 0.005, 
the solution becomes more singular. In order to fully resolve the spatial solution, we increase 
the spatial resolution to = 512. Following [23], we compute the time discretization error at 
time t as follows: 

eT(n; At) = \\u{T;At) - u{T; At/2)\\i2 . (209) 

For a vector field u(x) = (ni(x), U2(x)) defined on the Cartesian grid with Xi = ih, yj = jh, 
the discrete P norm is defined as follows 



\Mi^= \^('^^ix^^yj)+ul{xi,yj)jh^j . (210) 

Similarly, the discrete P norm for a vector field w(a) = {wi{a),W2{a)) defined on the interface 
r is defined below: 



= \y2[wliai)+w1{ai)) Aa] . (211) 




We compute the solution up to T = 1 and evaluate the convergence rate based on the numerical 
solution at T = 1 with different temporal resolutions. The results are shown in Fig. [7] and 
Table El As we can see, the convergence rate is approximately second order. 

Next we check the stability of the second order semi-implicit scheme. Fig [8] shows that 
the total energy for the second order explicit and second order semi-implicit schemes with 
different timesteps 0.002 and 0.02. We choose the same second order explicit scheme that was 
used in [13]. With At = 0.002, both the explicit and semi- implicit schemes are stable and 
they give nearly identical results. With At = 0.02, the explicit scheme becomes unstable, but 
the semi-implicit scheme is still stable. The stability restriction of the semi-implicit scheme is 
far less severe than the corresponding explicit scheme. Our numerical study also shows that 
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error of velocity field u 




-»^- |i=0.05 



Figure 7: Plot of P errors in time at time T = 1 for the second order semi-implicit scheme 
of the first kind. We choose 5;, = 1 and N = 256 in all computations except in the case of 
fi = 0.005 where N is increased to 512. The line at the bottom of each graph is a reference 
line which corresponds to the second order convergence rate. 





convergence rate of X 


convergence rate of u 


0.05 


2.11 


2.70 


0.01 


2.13 


2.10 


0.005* 


2.17 


1.96 



Table 3: Convergence rates for X and u fitted from the data shown in Fig [71 the case * is 
computed using a refined mesh 512 x 512. 




Figure 8: Total energy of the unsteady Stokes system for the second order semi-implicit scheme 
and the second order explicit scheme. = 128, Sf, = 1, fi = 0.01. Left: At = 0.002; Right: 
At = 0.02. 
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N 


/i = 0.05 


/i = 0.01 


fi = 0.005 


explicit 


semi-implicit 


explicit 


semi-implicit 


explicit 


semi-implicit 


64 


7 


0.8 


7 


0.8 


7 


0.8 


128 


37 


1.6 


38 


1.6 


38 


1.6 


256 


249 


4.4 


504 


4.6 


506 


4.5 


512 


3088 


24 


6182 


25 


6200 


25 



Table 4: Execution time for each computation in seconds. Here "explicit" stands for the second 
order explicit scheme and "semi-implicit" the second order semi-implicit scheme. 

the time step for the semi-implicit scheme is independent on the meshsize while the explicit 
scheme requires finer time step as the spatial mesh is refined. 

Finally, we compare the performance of our second order semi-implicit scheme with that 
of the second order explicit scheme. As before, in order to keep the accuracy with 5%, 
we take At = | for our semi-implicit schemes. For the explicit scheme, we take At = 
1/128,1/256,1/512,1/1024 which correspond to the spatial mesh sizes N = 64,128,256,512 
respectively, when fj, = 0.05. When ^ = 0.01 and 0.005, the time step is set to be At = 
1/128, 1/256, 1/1024, 1/2048. These time steps are the largest ones that we can take to keep 
the stability of the explicit scheme. We compute the solution up to T = 2. The result is shown 
in Table HI Again, we observe the same qualitative behavior as the first order schemes we 
reported earlier. 

8 Concluding Remarks 

In this paper, we have developed several efficient semi-implicit immersed boundary methods 
for solving the immersed boundary problem for the steady and unsteady Stokes equations. 
The immersed boundary method has emerged as one of the most useful numerical methods 
in computing fluid structure interaction, and has found numerous applications. But it also 
suffers from the severe time step stability limitation due to the stiffness of the elastic force. 
Guided by our stability analysis, we have developed several efficient semi- imp licit schemes 
which remove the stiffness of the immersed boundary method. We have demonstrated both 
analytically and computationally that our semi-implicit schemes have much better stability 
property than the explicit scheme. More importantly, unlike most existing implicit or semi- 
implicit schemes, our semi-implicit schemes can be implemented very efficiently. In fact, our 
semi-implicit scheme of the first kind has a computational cost that is essentially the same as 
that of an explicit scheme in each time step, but with a much better stability property. The 
saving in the computational cost is quite substantial. We have demonstrated this improved 
stability for a range of parameters and numerical resolutions. Our computational results show 
that the larger the spatial resolution is, the bigger the computational saving our semi-implicit 
schemes can offer. Thus the semi-implicit schemes we develop in this paper provide an effective 
alternative discretization to the explicit method. 

One of the essential steps in developing our semi-implicit schemes is to obtain an uncondi- 
tionally stable semi-implicit discretization of the immersed boundary problem. This provides 
us with a building block to construct our efficient semi-implicit schemes. There are two im- 
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portant observations in constructing the unconditionally stable semi-implicit discretization. 
The first one is that we need to preserve certain important solution structures at the discrete 
level. Specifically, we need to ensure that the discrete velocity field is divergence free, and 
the discrete spreading and interpolation operators are adjoint. Another essential step is to 
decouple the stiffness of the elastic force from the fluid flow in some appropriate way. This 
is difficult to achieve if we use the Cartesian coordinate. But it becomes easier if we use the 
arclength and tangent angle formulation to describe the dynamics of the immersed interface 
as was done in [9j. On the other hand, as we demonstrated in this paper, the unconditionally 
stable semi-implicit scheme is still very expensive to implement, and the saving over the expicit 
scheme is rather limited. 

Based on this unconditionally stable semi-implicit discretization, we have developed sev- 
eral efficient schemes for both the steady and the unsteady Stokes flows. By applying the 
Small Scale Decomposition to the unconditionally stable semi-implicit time discretization and 
further simplifying the leading order singular kernel, we obtain our semi-implicit scheme. The 
advantage of this semi-implicit scheme is that the leading order term can be expressed as 
a convolution operator, which can be evaluated explicitly using the Fourier transformation. 
This allows us to solve for the implicit solution explicitly in the spectral space, which offers 
substantial computational saving over the explicit scheme. 

It is a natural step to extend the the semi-implicit schemes developed for the unsteady 
Stokes equations to the Navier-Stokes equations. The discretization of the Navier-Stokes equa- 
tions shares many similar properties as the unsteady Stokes equations if we treat the convection 
term explicitly. We have performed a number of numerical experiments to test the stability 
and the robustness of our semi-implicit immersed boundary methods for the Navier-Stokes 
equations. The results are qualitatively similar to those for the unsteady Stokes equations 
which we have presented in this paper. These results will be reported in a subsequent paper. 
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A The semi-implicit scheme of the second kind 

In this appendix, we will derive the semi-implicit method of the second kind in more detail. 
As we mentioned before, the small scale decomposition only captures the leading order contri- 
bution from the high frequency components, which can not remove the stiffness induced by Sh 
and fi completely. The coefficients Sb and n can still affect the time stability through the low 
frequency components of the solution, which comes from the lower order term of the right hand 
side. In order to obtain a semi-implicit discretization with better stability property, we can 
incorporate the low frequency contribution from the second term in our implicit discretization. 
Eor the steady Stokes flow, this gives rise to the following decomposition: 

'-^-^ = -^|fc|5ri+.F(c/ln|a-a'|(.ri-l)Cda'), 
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CC/" -eij\n\a- a'\{sl - + 



4^' 



Ln+l 



At 



1 



n/in+1 



+ 



(212) 
(213) 



1 



where 7 = max 1 . By replacing the continuous derivative by the discrete derivative, 

"V SaJ 

and discretizing the continuous integral by the trapezoidal rule, we obtain our second semi- 
implicit scheme. We call this semi-implicit scheme the semi-implicit scheme of the second 
kind. Near equilibrium, we can show that the semi- implicit scheme of the second kind is 
unconditionally stable. Moreover, the stability property is independent of the meshsize, elastic 
coefficient 5^ and viscosity coefficient /u. Our numerical study also confirms this. The trade-off 
is that we need to solve a linear system at each time step to obtain the implicit solution at 

Similarly, in the case of unsteady Stokes flow, we can also include the second term of the 
right hand side in the leading order term to derive a scheme with better stability property. In 
this case, the leading order term becomes: 



T(s. 
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,jr f 1 yn+lgn+l 
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(215) 



where the derivative will be discretized by the spectral method and the integration will be 
discretized by the trapezoidal rule. We call the above scheme the semi-implicit scheme of the 
second kind for the unsteady Stokes flow. Near equilibrium stability analysis suggests that the 
semi-implicit scheme of the second kind is unconditionally stable. 



B Derivation of the semi-implicit scheme ( 11401) - ( fl43l) 



In this appendix, we will derive the semi-implicit scheme ()140p - (ll43p . We define the operator 
Q{sa'-, u", 0", X") : Sq ^ u by the following equations: 



u — u 



At 



-VhP + /uV|u + L^,„(F(s«, r)), 
V;,-L;,,„(F(s„,r)). 



(216) 

(217) 
(218) 
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Given Sa, we obtain u by solving above equations. From the definition of operator Q, we have 

u"+i =g(s^+i;u",r,X"), 
u*'"+i =g(s^;u",r,X"), 

where u*'""^^ is calculated from equations (|140p - (jl4ip . 
Then the equation of Sa can be rewritten as 

„n+l _ „n 

= DAa (^^,„(u"+i) • T^) - Z)A.r (l^,„(u"+i • n")) 
= Z)A.(L^,„(a(5r^u",r,X")).r") 

= 7^(s^+^u",r,X"). (219) 
Using the small scale decomposition, we get 

7^(sS+l; u", X") ~ T(sJJ+^). (220) 
By treating the leading order term implicitly, we obtain our semi-implicit method as follows: 

„n+l _ „n 

" " = ^(sr^) + (7^(5^;u^r,x-)-^(s-)) 



T(.r 1) + DAa {lIu {Q^sl; u", r , X")) • 

-i)A„r (g(.S; u", r , x")) • n") - t{sI) 



-r(.S) 

r(sS+i) + (DA„y*'"+^ - Z)Aae"t/*'"+' - T(s^)) . (221) 



This is exactly our semi- implicit scheme p34|) . 

In the steady Stokes case, we can define the operator ^(s^; u", 0", X") : — > u similarly: 

= -V^p + ^V|u + L;,,„(F(s„,r)), (222) 
VIp = Vh-Lh,n{F{sa,e^)). (223) 

In this case, we have 

g(s2;u",r,X") =u". 
Thus, the semi-implicit scheme becomes 

„n+l _ „n 

" " = T(5r') + (DAaV^ - DAa^f/- - T(5^)) . (224) 
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C The derivation of the Fourier transform of Ki 







In this appendix, we derive the the Fourier transform of Kq, which is given in (jl44p . By the 
definition of the Fourier transform, we have 

^ Ko{f3\a - a'\)f{a')da'^ 

Ko{(3\a - a'\)f{a')da'^ e-'^'^da 

Ko{P\a - a'\)f{a')e-'^''da'da 

Ko{^\a - a'|)/(a')e~^''^"""'^e-^'^'"'da'd(a - a') 

Ko{P\a - a'|)e-^^(°-°')(i(a - a')) f{a')e-'''"' da' 



1 r+oo 

J —OD \J — OO 

2 r+OD r+oo 



J— OO J — OO 
2 r+oo r+oo 



OO J — OO 



^ r+oo / r+oo 
VT 7-00 

= (225) 



OO 



where J-{f{a)){k) = / /(a)e* "'da is the Fourier transform. In the calculation above, we 

J— OO 

have used the expression of the Bessel function (p 376, [Ij) 



/ N /■+°° cos(te) , 

^o(^) = / -7=4dt' (226) 
Jo Vl + i 



and the identity below: 



OO JO \/T+l^ 



Ko{0\x\)e-'''''dx = I I ^^^^^^^^e-'^^dtdx 



1 /■+°° /■+~ cos(/3te) ,^ , 
' ' ^ -e '■''^dtdx 



^ J—oo J ~oo 



1 z'+oo c+oo piptx 

' ' -.e-'^''dtdx 



^ J —OO J—oo 
]^ r+oo r+oo ^i{l3t—k)x 



rdtdx 



^ J—oo J—oo 

= .r'^^lZ^dt= , - ■ (227) 

This proves ()144p . 
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